home *** CD-ROM | disk | FTP | other *** search
/ Aminet 2 / Aminet AMIGA CDROM (1994)(Walnut Creek)[Feb 1994][W.O. 44790-1].iso / Aminet / util / misc / Fudgit233.lha / Source / examples / advance.ft next >
Encoding:
Text File  |  1993-12-14  |  4.1 KB  |  182 lines
  1. # This file is a real calculation.
  2. # It builds a Wulff shape for a 2-dimensional Ising model
  3. # using Onsager solution of the surface tension.
  4. # The final plot is the shape of a droplet together with
  5. # the polar plot of the surface tension
  6. # It is square at low temperature and circular at higher T.
  7. # T are in units of Tc.
  8. # Calculations using expressions from:
  9. # Avron et al., J. Phys. A 15, L81 (1982).
  10. # Never mind.....
  11. #
  12. # Plot temperature dependence of epsilon
  13. # fitting preamble
  14. pmode
  15.     # set term uniterm
  16.     # set term X11
  17.     # set term postscript
  18.     # set output 'fig.ps'
  19.     # set size 3.5/5.0, 3.5/3.5
  20.     set grid
  21.     set nokey
  22.     set data style line
  23.     # set xrange [0:1.1]
  24.     set auto
  25.     set format '% .1lf'
  26.     set polar
  27. fmode
  28. # define variables
  29. # Critical temperature of the Ising model
  30. let tc = -2.0/ln(sqrt(2)-1)
  31. set data 50
  32.  
  33. # A macro to take rough derivative 
  34. # Syntax deriv X Y
  35. # creates vector DY
  36. macro deriv 2
  37.     cmode
  38.         for (i=2;i<data;i++) {
  39.             D$2[i] = ($2[i+1] - $2[i-1])/($1[i+1] - $1[i-1])
  40.         }
  41.         D$2[1] = D$2[2]
  42.         D$2[data] = D$2[data-1]
  43.     fmode
  44. stop
  45.  
  46. # build angular vectors from 0 to pi/2
  47. cmode
  48.     n=1
  49.     tmp = data + 1
  50.     PHI = n++/tmp
  51.     PHI *= pi/2
  52.     THETA = PHI
  53.     CF = cos(PHI)
  54.     CF2 = CF^2
  55.     C2F = cos(2 * PHI)
  56.     C2F2 = C2F^2
  57.     SF = sin(PHI)
  58.     SF2 = SF^2
  59.     S2F = sin(2 * PHI)
  60.     S2F2 = S2F^2
  61. fmode
  62.  
  63. # polar vector WULF to be minimized
  64. # with the following cmode function
  65. cmode
  66.     func minim(x, y) {
  67.         if (x < y) {
  68.             return(x)
  69.         }
  70.         return(y)
  71.     }
  72. fmode
  73.             
  74. # define macro for temperature dependence
  75. macro phimake 2
  76.     echo Doing T = $1
  77.     cmode
  78.     # build scalars
  79.         tt = $1
  80.         t = tt*tc
  81.         x = exp(-2/t)
  82.         x2 = x * x
  83.         dx = -2 * x
  84.         a = (1 + x2)^2
  85.         a2 = a^2
  86.         da = -8 * x2 * (1 + x2)
  87.         p = 2 * x * (1 - x2)
  88.         p2 = p^2
  89.         dp = -4 * x * (1 - 3 * x2)
  90.     # build vectors
  91.         B = sqrt(0.25 * a2 * S2F2 + p2 * C2F2)
  92.         DB = (0.25 * a * da * S2F2 + p * dp * C2F2)/B
  93.         NUM1 = (a2 * SF2 + p2 * C2F)
  94.         DNM1 = (a * SF2 + B)
  95.         AL1 = acosh(NUM1/(DNM1 * p))
  96.         NUM2 = (a2 * CF2 - p2 * C2F)
  97.         DNM2 = (a * CF2 + B)
  98.         AL2 = acosh(NUM2/(DNM2 * p))
  99.         DNUM1 = 2 * (a * da * SF2 + p * dp * C2F)
  100.         DDNM1 = (da * SF2 + DB)
  101.         DNUM2 = 2 * (a * da * CF2 - p * dp * C2F)
  102.         DDNM2 = (da * CF2 + DB)
  103.         TMP1 = (DNUM1 * DNM1 - NUM1 * DDNM1)/(DNM1 * DNM1)
  104.         TMP2 = (DNUM2 * DNM2 - NUM2 * DDNM2)/(DNM2 * DNM2)
  105.         DAL1 = (TMP1 - dp * cosh(AL1))/(p * sinh(AL1))
  106.         DAL2 = (TMP2 - dp * cosh(AL2))/(p * sinh(AL2))
  107.         S = (AL1 * SF + AL2 * CF) * t
  108.         E = DAL1 * SF + DAL2 * CF
  109.     # build Wulff vector
  110.         WULF = S
  111.         for (i=1;i<=data;i++) {
  112.             for (j=1;j<=data;j++) {
  113.                 TMP[j] = S[j]/cos(PHI[i] - THETA[j])
  114.                 WULF[i] = minim(WULF[i], TMP[j])
  115.             }
  116.         }
  117.     fmode
  118.     # take derivative -> creates vector DWULF
  119.     deriv PHI WULF
  120.     # recalculate E(BETA) to average
  121.     cmode
  122.         BETA = PHI - atan(DWULF/WULF)
  123.         CA = cos(BETA)
  124.         CA2 = CA^2
  125.         C2A = cos(2 * BETA)
  126.         C2A2 = C2A^2
  127.         SA = sin(BETA)
  128.         SA2 = SA^2
  129.         S2A = sin(2 * BETA)
  130.         S2A2 = S2A^2
  131.         B = sqrt(0.25 * a2 * S2A2 + p2 * C2A2)
  132.         DB = (0.25 * a * da * S2A2 + p * dp * C2A2)/B
  133.         NUM1 = (a2 * SA2 + p2 * C2A)
  134.         DNM1 = (a * SA2 + B)
  135.         AL1 = acosh(NUM1/(DNM1 * p))
  136.         NUM2 = (a2 * CA2 - p2 * C2A)
  137.         DNM2 = (a * CA2 + B)
  138.         AL2 = acosh(NUM2/(DNM2 * p))
  139.         DNUM1 = 2 * (a * da * SA2 + p * dp * C2A)
  140.         DDNM1 = (da * SA2 + DB)
  141.         DNUM2 = 2 * (a * da * CA2 - p * dp * C2A)
  142.         DDNM2 = (da * CA2 + DB)
  143.         TMP1 = (DNUM1 * DNM1 - NUM1 * DDNM1)/(DNM1 * DNM1)
  144.         TMP2 = (DNUM2 * DNM2 - NUM2 * DDNM2)/(DNM2 * DNM2)
  145.         DAL1 = (TMP1 - dp * cosh(AL1))/(p * sinh(AL1))
  146.         DAL2 = (TMP2 - dp * cosh(AL2))/(p * sinh(AL2))
  147.         EOFB = DAL1 * SA + DAL2 * CA
  148.         tote = totl = 0
  149.         for (i=1;i<=data;i++) {
  150.             tote += WULF[i] * EOFB[i]
  151.             totl += WULF[i]
  152.         }
  153.         tote /= totl
  154.     fmode
  155.     save vec PHI E S WULF BETA $Tmp.$2
  156.     append var tt tote /tmp/ener.wulf
  157. stop
  158.  
  159. # make a few
  160. ! rm -f /tmp/ener.wulf
  161. phimake 0.01   0
  162. let temperature = 0.1
  163. let plotnum = 1
  164. while (temperature <= 0.9)
  165.     phimake $temperature    $plotnum
  166.     let plotnum++
  167.     let temperature+=0.1
  168. end
  169. phimake 0.95  10
  170. !sync
  171.  
  172. pmode
  173. #    plot '$Tmp.0', '$Tmp.1', '$Tmp.2', '$Tmp.3' , '$Tmp.4', '$Tmp.5'
  174. #    plot '$Tmp.0' us 1:3, '$Tmp.1' us 1:3, '$Tmp.2' us 1:3, \
  175. #    '$Tmp.3' us 1:3, '$Tmp.4' us 1:3, \
  176. #    '$Tmp.0', '$Tmp.1', '$Tmp.2', '$Tmp.3' , '$Tmp.4'
  177.     plot '$Tmp.0' us 1:3, '$Tmp.0' us 1:4,\
  178.     '$Tmp.4' us 1:3, '$Tmp.4' us 1:4,\
  179.     '$Tmp.7' us 1:3, '$Tmp.7' us 1:4
  180. fmode
  181.  
  182.